RNA can carry out many biological functions, such as translation, transcription, conservation of genetic information, and regulation of gene expression.^[1–4] To perform these biological functions, RNA adopts a variety of complex three-dimensional structures.^[5–9] The functions of many RNAs are often kinetically controlled.^[10–16] For instance, self-induced riboswitches can regulate gene expression by limiting the folding of functional structures to certain time windows.^[15,16] The hok/sok system of plasmid R1^[16] regulates the plasmid maintenance through mRNA conformation rearrangements into different functional forms. A nucleotide base pair is the basic unit of nucleic acids molecular structures. In the process of structural rearrangement, closing and opening of base pairs are fundamental events. So studying the thermodynamic and kinetic properties of opening and closing of a single base pair is crucial for further exploring and predicting the biological functions of RNA molecules.^[17–20]

By assuming the nearest-neighbor interaction, the thermodynamic parameters of a base stacking have been experimentally determined by measuring the melting curves of a few duplexes with several base pairs through ultraviolet absorption.^[21] The accuracy of the thermodynamic parameters of a base pair depends on the assemble of the duplexes and a great correction has been made since its first version.^[21] The exchange of imino protons^[22] and NMR^[23–26] have been used to study the kinetic properties of single nucleotide base pair. However, due to the randomness of base pair folding, the complexity of overcoming the energy barrier,^[27–29] as well as the limitations and effects of the experimental instrument,^[30–32] it is very difficult to directly observe the detailed kinetic process and conformational changes of the base pair opening and closing in experiments. Molecular dynamic simulation, which can make up the experimental deficiencies, has been widely used to study the bio-molecules.^[33–45] Recently, molecular dynamic simulation has been used to study the kinetic process of the DNA and RNA base pair opening and closing. Not only the transition states and pathways between open and closed states of the nucleotide base pair were identified, but also the thermodynamic parameters and the dynamic rates of a base pair opening and closing were obtained, which are in good agreement with those of the experiments. Here we will provide a brief overview of the studies about the thermodynamic and opening-closing kinetic properties of a nucleotide base pair.

Using four parameters, namely the distance d₁₆ and interaction energy E₁₆ between the terminal base pair G6–C1, as well as the distance d₁₂ and interaction energy E₁₂ between the adjacent base G2–C1, Hagan et al.^[39] have studied the kinetic pathways for DNA terminal single base pair G6–C1 opening and closing through molecular dynamic simulation. The closing and opening of the base pair were defined by the four parameters. When the two interaction energy parameters were less than a certain value ( and ), the corresponding structures could be regarded as the closed state, and when the two distance parameters were within a certain range ( and ), the corresponding structures could be regarded as the open state. Two qualitatively different pathways of the terminal base pair G6–C1 from the closed state to open state were identified through this method. One path was that the flipping base first opened its intramolecular hydrogen bonds of complementary base pair before it was unstacked, and the other path was that both sets of base stacking interaction between adjacent bases were ruptured simultaneously. The two kinetic pathways rejoined at the transition states, where the magnitudes of the interaction energies of E₁₂ and E₁₆ were comparable with the thermal energy.

Using atomistic pulling simulations and reweighting scheme, Colizzi et al.^[45] have studied the forming and melting of four different terminal complementary base pairs G–C in RNA duplexes. The free-energy profile was reconstructed as a function of the number of water molecules surrounding the unbinding base, which could be used to distinguish the configuration from the other and to properly collect the corresponding weights. By fitting the free-energy profiles with the combination of two quadratic functions, the free energy difference between the open and closed states was calculated. According to a systematic step-by-step approach, two different opening paths from the closed state to open state were identified. Each path could be divided into two steps. First, the complementary base pair was partially opened by the unpaired base of the 5-terminal or 3-terminal; second, the remaining dangling base at the other terminal was then unstacked, and the terminal base pair was completely opened. The free energy changes from the closed state to the relatively stable intermediates and from the intermediates to the open state in each path were attained.

Recently, Xu et al.^[42] applied a mixed method, which integrated molecular dynamic simulation, kinetic Monte Carlo simulation and master equation methods, to study the kinetic mechanism of RNA single base pair formation. The molecular dynamic simulation method was used to get the trajectory including many structures. All the conformations were assigned to 50 clusters according to the structural closeness by the RMSD (Root Mean Square Deviation), as well as three order parameters: the distance d₁₂ between the geometric centers of all the heavy atoms in the G1 and G2 bases and the nonbonded interaction energies E₁₂ and E₁₆ between the sequentially neighboring nucleotides G1 and G2 and between the pairing nucleotides G1 and C6, respectively. Based on the behavior of the time-dependent population for each cluster, the 50 clusters were classified as four states: unfolded state, intermediate state, trapped state and folded state. At last, the overall folding kinetics could be described by the four-state kinetics scheme.

Since it usually takes a long time to open or close the RNA base pair at room temperature,^[46] the simulation near its melting temperature of the base pair can quickly make the base pair reach the two-state equilibration between open and closed states. The melting temperature of the AU base pair can be estimated as: , where and are the enthalpy and entropy changes of the terminal base-stacking from the nearest-neighbor model.^[47] The thermodynamic and kinetic properties of RNA terminal AU^[48] base pair were attained through molecular dynamic simulation at different temperatures: 390 K, 400 K, 410 K, 420 K, and 430 K, nearing the melting temperature of the AU base pair. The initial structure of a 12-nucleotides oligonucleotide -(AAGCGCAAGCUU)- , which includes a tetraloop and four base pairs, was obtained from its crystal structure (PDB ID:1ZIH). To specifically study the thermodynamic parameters and kinetic mechanism of the opening and closing of the terminal AU base pair, all other nucleotides except the two terminal nucleotides -A and -U of the RNA hairpin were fixed by using an additional force.

According to the simulated trajectories, the time-dependent RMSD of the two terminal nucleotides -A and -U relative to the initial structure, where i represents the i-th atom, the total number of atoms, and represent the coordinates of the atom i at time zero and t, respectively, and the torsional dihedral angles ζ (C ,^[49] where j represents the j-th nucleotide in the polynucleotide chain) on the backbone of the terminal base pair AU were obtained. These values clearly showed a two-state distribution, as shown in Fig. 1. Therefore, according to the distribution of RMSD and torsional angles, we divided the simulation structures into closed state (cs), open state (os), and transition state (ts). When the RMSD value was less than 2.0 Å, at which the corresponding torsional angle ζ was always between −50° and −100°, the conformation was treated as the closed state. On the contrary, when the RMSD value was larger than 2.0 Å and the corresponding torsional angle ζ was always between 50° and 100°, the conformation was considered to be the open state (see Figs. 2(a) and 2(b)). When RMSD values were greater than 2.0 Å, but the torsional angle ζ were still between −50° and −100°, then such conformations were regarded as the transition state (see Figs. 2(a) and 2(b)). The trajectories of transition paths showed that all the transitions from the closed states to open states or from the open states to closed states go through the transition states. Namely, in the process base pair from the closed states to open states, the bases of terminal nucleotides first flipped out from the closed states with the increased RMSD, and then the torsional angles of backbone changed from the range of (−50°)–(−100°) to the range of (50°)–(100°), however, in the process of base pair from the open states to closed states, the torsional angles of backbone first changed from the range of (50°)–(100°) to the range of (−50°)–(−100°), and then the flipped bases of terminal nucleotides returned to the closed states with the reduced RMSD. According to the location of the transition state appeared, the transition state could be classified into two types: ctc and oto, where ctc referred to the transition state that transited from the closed state and then back to the closed state, and oto referred to the transition state that transited from the open state and then back to the open state (see Fig. 2(c)).

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The RMSD values and torsional angles ζ during the whole simulation times ( ) at 400 K.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) (a) The RMSD values among the simulation times of 412 ns∼463 ns at 410 K. (b) The torsion angles ζ for the simulation times of 412 ns∼463 ns at 410 K. (c) The RMSD and torsion angles ζ near the transition states.

As compared to the open state and closed state, the probabilities of the transition state were very small. According to the temperature dependence population distribution of the closed state and open state, the enthalpy change and the entropy change of the terminal AU base pair opening and closing could be attained as (see Fig. 3), where and are the occupied population of the open states and closed states, respectively; is the Boltzmann constant; T is the absolute temperature. The results showed that thermodynamic parameters were in great agreement with the nearest-neighbor model.^[47]

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The temperature dependence of . Symbol is the simulation results and line is the linear fit.

The average lifetime of each state could be calculated as , where is the average lifetime, N is the corresponding total number of times the molecules reside in the closed state or open state or transition state, and τ_i is the i-th lifetime of the conformation in the corresponding state. The average lifetime of the closed state and open state showed different temperature-dependent behaviors. The average lifetime of the closed state showed the strong temperature-dependent behavior, while that of the open state demonstrated the weak temperature dependence behavior, as shown in Fig. 4(a). Compared to the open state and closed states, the average lifetimes of transition states ctc and oto were very short and also exhibited the weak temperature dependence behavior, as shown in Table 1.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (a) Temperature dependence of the average lifetime of the closed state (hollow triangle) and the open state (hollow circle). Symbol is the simulated results. Lines are fitted results with Eqs. (1) and (2). (b) Temperature dependence of the ratios of (hollow triangle) and (hollow circle). Symbol is the simulated results. Lines are fitted results with Eqs. (1)–(4).

Table 1.

Table 1.

Table 1. The average lifetime τave (in unit ns), the occupied probability p, and the total number N of occurrences of conformations at closed, open, and transition states at different temperatures (in unit K) during the 3000-ns simulation time. . Temperature/K Open state (os) Closed state (cs) Transition state (oto) Transition state (ctc) T No po Nc pc Noto poto Nctc pctc 390 24.207 58 0.468 27.051 59 0.532 0.295 79 0.0078 0.309 212 0.0218 400 23.612 67 0.527 20.853 68 0.473 0.288 104 0.0100 0.292 304 0.0296 410 23.187 75 0.580 16.813 75 0.420 0.284 153 0.0145 0.284 391 0.0370 420 22.783 83 0.630 13.202 84 0.370 0.280 185 0.0173 0.282 335 0.0315 430 22.451 91 0.681 10.402 92 0.319 0.278 271 0.0251 0.281 363 0.0340

Table 1. The average lifetime τave (in unit ns), the occupied probability p, and the total number N of occurrences of conformations at closed, open, and transition states at different temperatures (in unit K) during the 3000-ns simulation time. .

For the opening-closing two-state transition, the closing rate and the opening rate of the base pair could be calculated from the average lifetime of the open and closed states as , , where and are the average lifetime of the open and closed states, respectively. The transition rates from the transition state to the closed state ( ) and to the open state ( ) could be calculated according to the average lifetime of the transition state ctc ( ) and oto ( ) based on the following formula, and . The transition path time from the closed state to the open state , and from the open state to the closed state , were both weakly temperature-dependent behavior.

According to the transition-state theory,^[50–54] the average residence time in the closed state ( ) and open state ( ), the transition path time and could be expressed by the following equations 1 2 3 4 where and are free energy barrier heights of the closed state and open state; ( ) is the diffusion coefficient at the free energy barrier top; , , and are the curvatures of the free energy surface at the barrier top, the closed state, and the open state, respectively. According to Eqs. (1)–(4), we found that showed strong temperature-dependent behavior, while the showed temperature-independent behavior, as shown in Fig. 4(b).

Assuming that was temperature-independent as in protein folding,^[54] the free energy barrier of closing the base pair from the open state should depend on temperature as . By fitting the above two curves, we found that the free energy barrier from the closed state to the open state corresponded to the enthalpy change, while the free energy barrier of the reverse transition corresponded to the entropy change. The temperature dependence of the residence time in the open state came from the prefactor of Eq. (2), which demonstrated that the simple Einstein relationship , where η is the viscosity of the solvent, was not applicable to base folding. Assuming that as for protein folding, in which is the local mean-squared fluctuation in energy and is a measurement of the underlying landscape roughness, was found to be . These results suggested that a one-dimensional free energy surface was sufficient to accurately describe the dynamics of base pair opening and closing, and the dynamics were Brownian.

In conclusion, by simulating the opening-closing switch transition of RNA single base pair, we could not only qualitatively obtain the thermodynamic and kinetic properties of the base pair, such as the kinetic transition pathways, the transition states and the kinetic mechanisms of the base pair opening and closing, in which a one-dimensional free energy surface was sufficient to accurately describe the dynamics of RNA single base pair opening and closing, and the dynamics were Brownian. But also the quantitative thermodynamic parameters, such as the enthalpy and entropy changes, and the kinetic parameters, such as the transition rates and the transition path time, were in good agreements with the experimental results upon the base pair opening and closing. Although molecular dynamics simulation has its limits, such as the force constants and simulation times,^[55] it would be a powerful tool to further study the thermodynamic and kinetic properties of base pair and base stacking interaction,^[56,57] the next nearest-neighbor and the nearest-neighbor effects,^[58] as well as pseudoknot,^[59] and the ion effects,^[60] which are still hard to experimentally determine, but are important to understand the biological relevance of RNA structures.